# REPLICATION INFO ####
# prepared for replication April 25, 2023
# R version 4.2.2 (2022-10-31) -- "Innocent and Trusting"
# macOS Ventura 13.0.1

# packages ####
library(tidyverse)
library(gridExtra) 

# read in data ##### 
# (skip first row which provides question labels)
dat <- (read.csv("ethics_apsa_survey.csv")[-1,])

names(dat)

# EVALUATION RANKINGS ####

## convert to numeric
## variables rank_ethics (7th) through harm_science_labexp (33rd)
dat_num <- as.data.frame(dat[,7:33] <- 
                             lapply(dat[,7:33,drop=FALSE],as.numeric))

dat_ranks <- dat_num %>% 
  rename(
    Ethics=rank_ethics,
    Writing=rank_writingquality,
    Design=rank_designquality,
    Analysis=rank_analysisquality,
    Theory=rank_theory,
    Results=rank_results
  )

## reshape 
dat_ranks <- dat_ranks %>% pivot_longer(Ethics:Results)
# order "name" (item) by mean "value" (rank)
(avgrank <-aggregate(value ~ name, dat_ranks, mean))
dat_ranks$name <-factor(dat_ranks$name, levels=avgrank[order(avgrank$value), "name"])

## plot evaluation rankings (Figure 1) ####
rank_plot <- ggplot(dat_ranks, aes(x=name,y=value)) +
  stat_summary(fun.data= mean_cl_normal, fun.args = list(conf.int = 0.95)) + theme_bw() + xlab("") + 
  ylab("Less Important                         More Important") +
  coord_cartesian(ylim=c(6, 1)) + 
  theme(axis.text=element_text(size=14),axis.title=element_text(size=14))+
  scale_y_continuous(breaks = c(1:6))

rank_plot

## export plot ####
pdf("rankings.pdf", width=6, height=5)
rank_plot
dev.off()

# PERCEIVED POTENTIAL HARM ####

## harm to subjects ####
dat_harm1 <- dat_num %>% 
  rename(
    Audit=harm_subjects_audits,
    Ethnography=harm_subjects_ethnog,
    Field.Exp=harm_subjects_fieldexp,
    Survey.Exp=harm_subjects_surveyexp,
    Lab.Exp=harm_subjects_labexp,
    Interviews=harm_subjects_interviews,
    Observational=harm_subjects_observ
  )

# reshape
dat_harm1 <- dat_harm1 %>% pivot_longer(Interviews:Lab.Exp)

# order "name" (item) by mean "value" (rank)
(avgharm1 <-aggregate(value ~ name, dat_harm1, mean))
dat_harm1$name <-factor(dat_harm1$name, levels=avgharm1[order(avgharm1$value), "name"])

### plot harm to subjects (left plot Figure 2) ####
harm1_plot <- ggplot(dat_harm1, aes(x=name,y=value)) +
  stat_summary(fun.data= mean_cl_normal, fun.args = list(conf.int = 0.95)) + theme_bw() + xlab("") +
  ylab("Less Harmful                         More Harmful") +
  coord_cartesian(ylim=c(6, 1)) + 
  theme(axis.text.x = element_text(angle = 35, hjust = 1, size=11),axis.title=element_text(size=12),plot.title = element_text(size = 14)) +
  scale_y_continuous(breaks=c(1:7)) +
    ggtitle("Research Subjects")

harm1_plot

## harm to society ####
dat_harm2 <- dat_num %>% 
  rename(
    Audit=harm_society_audits,
    Ethnography=harm_society_ethnog,
    Field.Exp=harm_society_fieldexp,
    Survey.Exp=harm_society_surveyexp,
    Lab.Exp=harm_society_labexp,
    Interviews=harm_society_interviews,
    Observational=harm_society_observ
  )

# reshape
names(dat_harm2)
dat_harm2 <- dat_harm2 %>% pivot_longer(Interviews:Lab.Exp)

# keep order from first plot
(avgharm2 <-aggregate(value ~ name, dat_harm2, mean))
dat_harm2$name <-factor(dat_harm2$name, levels=avgharm1[order(avgharm1$value), "name"])

### plot harm to society (middle plot Figure 2) ####
harm2_plot <- ggplot(dat_harm2, aes(x=name,y=value)) +
  stat_summary(fun.data= mean_cl_normal, fun.args = list(conf.int = 0.95)) + theme_bw() + xlab("") +ylab("") + 
  coord_cartesian(ylim=c(6, 1)) + 
  theme(axis.text.x = element_text(angle = 35, hjust = 1, size=11),axis.title=element_text(size=12),plot.title = element_text(size = 14)) +
  scale_y_continuous(breaks = c(1:7)) +
  ggtitle("Society")

harm2_plot

## harm to scientic community ####
dat_harm3 <- dat_num %>% 
  rename(
    Audit=harm_science_audits,
    Ethnography=harm_science_ethnog,
    Field.Exp=harm_science_fieldexp,
    Survey.Exp=harm_science_surveyexp,
    Lab.Exp=harm_science_labexp,
    Interviews=harm_science_interviews,
    Observational=harm_science_observ
  )

# reshape
names(dat_harm3)
dat_harm3 <- dat_harm3 %>% pivot_longer(Interviews:Lab.Exp)

# keep order from first plot
(avgharm3 <-aggregate(value ~ name, dat_harm3, mean))
dat_harm3$name <-factor(dat_harm3$name, levels=avgharm1[order(avgharm1$value), "name"])

### plot harm to scientific community (right plot Figure 2) ####
harm3_plot <- ggplot(dat_harm2, aes(x=name,y=value)) +
  stat_summary(fun.data= mean_cl_normal, fun.args = list(conf.int = 0.95)) + theme_bw() + xlab("") +ylab("") + 
  coord_cartesian(ylim=c(6, 1)) + 
  theme(axis.text.x = element_text(angle = 35, hjust = 1, size=11),axis.title=element_text(size=12),plot.title = element_text(size = 14)) +
  scale_y_continuous(breaks = c(1:7)) +
  ggtitle("Scientific Community")

harm3_plot

## export combined plots ####
pdf("rankings_harmful.pdf", width=12, height=4.5)
grid.arrange(harm1_plot, harm2_plot, harm3_plot, ncol=3)
dev.off()

# DESCRIPTIVE STATS OF SAMPLE (Appendix) ####
library(sjPlot)

a <- as.data.frame(prop.table(table(dat$position2)))
b <- as.data.frame(prop.table(table(dat$subfield)))
c <- as.data.frame(prop.table(table(dat$gender)))
demos <- rbind(a, b, c)
demos <- demos[demos$Var1 != "", ]
demos$Freq <- demos$Freq*100
colnames(demos) <- c(" ", "Percent of Sample (%)")
tab_df(demos, file="survey_demos.doc")



